A solvable senescence model showing a mortality plateau 



O 

o 

(N 

> 
O 

(N 



J. B. Coe 1 , Y. Mao 1 , M. E. Cates 2 
1 Cavendish Laboratory, Madingley Road, Cambridge, CBS OHE, UK 
2 Dept. of Physics and Astronomy, University of Edinburgh, 
King's Buildings, May field Road, Edinburgh, EH9 3JZ, UK 
(Dated: February 1, 2008) 

We present some analytic results for the steady states of the Penna model of senescence, gen- 
eralised to allow genetically identical individuals to die at different ages via an arbitrary survival 
function. Modelling this with a Fermi function (of modest width) we obtain a clear mortality plateau 
late in life: something that has so far eluded explanation within such mutation accumulation mod- 
els. This suggests that factors causing variable mortality within genetically identical subpopulations, 
which include environmental effects, may be essential to understanding the mortality plateau seen 
in many real species. 

PACS numbers: 87.10.+e, 87.23.-n 
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A common feature of the life-history of multicellular 
organisms is the progressive decline in the various physio- 
logical processes once the reproductive phase is complete. 
This process of aging has attracted considerable atten- 
tion H, §, §, |, §, particularly after work from Human 
Genome project identified specific 'clock genes' Q which 
regulate human aging. In addition to these underlying 
genes the environment can play a part. For example, 
experimental studies |t], || have shown that organisms 
subject to a reduction in caloric intake without essen- 
tial nutrient deficiency display an extension of maximum 
life-expectancy. How, then, do the effects of evolution 
and/or the environment 'pick out' a particular profile of 
senescence? 

Two major theories of aging prevail, the 'mutation 
accumulation' where mutations affecting old ages accu- 
mulate thanks to a weaker reproductive selection, and 
the 'antagonistic pleiotropy' where a gene, beneficial in 
youth, is deleterious in old age Q . Traditional mutation 
accumulation theories lead to an exponentially increas- 
ing mortality with age, the so-called Gompertz Law |Q. 
However, it has been experimentally observed that the 
rate of mortality, whilst obeying the Gompertz law up to 
an intermediate age, often shows an unexpected drop at 
an advanced age. This drop, valid for many species from 
human to mcdflics, gives rise to a 'mortality plateau', 
somewhat at odds with theory Jj| |l^, 12 . As a re- 
sult, 'antagonistic pleiotropy' has been suggested as an 
essential part of aging j|, |lj . 

In 1995, with computer simulation in mind, Penna 
proposed a mutation accumulation model of senescence, 
which quickly became widespread | [L4| . Gompertz be- 
haviour has been reported for the standard Penna model 
fL5[ , and variations on the model were proposed to ac- 
count for demographic features |H| [l^, jig] , as well as 
catastrophic mortality | 20[|. H owever, a major short- 
coming of the Penna modei pifl was examined only re- 
cently: according to the model, all members of a genet- 
ically identical population must die exactly at the same 
age, which is obviously untrue. Recent attempts in re- 
moving this determinism within the Penna model gave 



rise to a gently decelerating old age mortality, without 
reproducing a satisfactory 'mortality plateau' p2| . 

Here we present a general formulation and an exact 
solution to the Penna model, which together with a sim- 
ple Fermi survival function, gives rise to a convincing 
mortality plateau. The introduction of a gene-dependent 
survival curve (replacing a single parameter, the deter- 
ministic age of death) enlarges the parameter space for 
simulations considerably. This makes analytical results 
all the more valuable, especially those that do not depend 
on the survival function chosen. 

Below, we first consider a particular version of the 
Penna model (by setting the mortality threshold T = 1), 
and present for it an exact analytic solution. Further- 
more, this solution is robust against changes in the sur- 
vival function. In the case of a Fermi function for sur- 
vival, we show the presence of an extended 'mortality 
plateau'. The generalisation of the analytic solution to 
T > 1 is then given. 

In the original Penna model E3, each individual car- 
ries a string of binary numbers, which stays fixed for that 
individual's lifetime. A 1 for the i-th bit represents the 
effect of a harmful gene which causes the individual to be 
struck by a heritable disease upon reaching the age of i. 
Thus, as an individual ages, its bit-string is sequentially 
examined and diseases are accumulated. The individual 
expires upon encountering the T-th deleterious bit, where 
T is a preset mortality threshold. Qualitatively, T repre- 
sents the number of heritable diseases required to make 
an individual nonviable; within the original model, if the 
T-th 1 occurs at the /-th bit, the individual will survive 
precisely to age I and no more. As long as it survives, 
the individual has a probability (rate) b to give birth. 
The offspring inherits the same bit-string from the par- 
ent except for some mutations occuring at a small rate 
(3. In simulations, the birth rate b is often regulated by 
a 'Verhulst factor' to prevent either population explosion 
or decay [Q. A different Verhulst factor, reducing popula- 
tions of all ages at every time step (but with a fixed birth 
rate), has also been used, for which the Penna model can 
be formulated analytically into an eigenvalue equation 
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and solved numerically p3fl . 

In the T = 1 limit, the population can be described 
by a probability function n(l) of the first deleterious bit 
occuring at position I. Equivalently n(l) is simply the 
number of such Z-type individuals at any particular time. 
Then the evolution of n(l) is given by 



dn(l) 
dt 



I' =1+1 



n(l') (1) 



where the first term corresponds to the mutation-free re- 
production; the second term gives the mortality, where 
l/l is the average rate of mortality for an /-type (which 
assumes the population changes slowly in the lifetime of 
an individual); the third term represent the switching 
from I' to I due to a single mutation p^| . Note that 
there is no noise in this equation. Hence it describes 
the thermodynamic limit of a large population in which 
deterministic dynamics is recovered. Note also that in 
a stationary state of the original Penna model, Z-types 
are distributed uniformly over all ages up to / at which 
they promptly die, so the mortality rate is exactly l/l. 
However, our formulation is more general in that Z-types 
may have an arbitrary survival function provided their 
mortality rate averages to 1/7. (This is the same as de- 
manding the average life expectancy for an Z-type to be 
I, which we may take as the definition of I.) As in the 
original Penna model, only bad mutations (from to 1) 
have been included. 

For a stationary state, the constant introduction of bad 
mutations is balanced by the longer reproductive lifetime 
of healthier individuals. Setting dn/dt = 0, Eq.(Q) can 
be solved exactly to give the recursion relation: 



n(l + 1) 1+ 1 



s# - bl 



n(l) 



I ePV+V -b(l + l)e-P 



(2) 



Before going further, we need to examine the inter- 
dependence of the mutation and birth rates (3 and b. Ob- 
viously a very large b coupled with a negligible (3 would 
result in an exponential population growth. A station- 
ary state requires a sepcific combination of (3 and b. An 
extremely long-lived individual (large I) producing many 
offspring in total (proportional to I), but mutation re- 
duces the probability of it reproducing itself accurately 
by an exponential factor (e^ 131 ). This means that a very 
large I cannot maintain itself, and must rely instead on 
mutated reproduction of even fitter individuals. Such a 
cascade is not possible in a finite population which must, 
at any time, contain a maximum I = l ma x] the above 
argument shows that, if too large initially, this will de- 
crease with time until the l ma x subpopulation is indeed 
self sustaining. The evolution equation Eq.(|l]) for this 
subpopulation then reads: 



dn(l rl 



dt 



dOmax^) Tl(lmax) 



where the growth rate function g(l) = be @ l — l/l, and 
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FIG. 1: Lifespan distribution n(Z) for l max = 20 (+), 30 (x), 
compared with simulations (boxes). Simulation size 10 7 , av- 
eraged over 10 runs. 



steady state demands g(l r . 



,-pi max 



1/lr 



0, namely: 



= 0. 



(4) 



Next we argue that for a steady state in the thermody- 
namic limit, we should have: 



I On 



1)<0; g(l ma x + 1) < 0. 



(5) 



Without the first condition, the subpopulation with I = 
lmax — 1 would grow thanks to the extra proliferation of 
lmax — 1 as a result of a mutation to l max ■ The second con- 
dition ensures the stability of l max + 1 upon approaching 
steady state. In the continuum limit, this is equivalent 
to demanding that the derivative of g{l) should be zero. 

These conditions are well confirmed by our numerical 
simulations of the Penna model. However, if n(l max + 1) 
falls to zero discontinously significantly before the steady 
state is reached, then the second condition in Eq.(||) 
might not be strictly required. This is unlikely in the 
thermodynamic limit and we ignore such contingencies 
from now on. It would also happen if a maximum bit- 
string length less than l max were imposed in a simulation; 
see below. 

Equations (f|,||) together lead to: 
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1 



(6) 



which implies a unique integer value for l ma x given any f3. 
Whenever 1/(3 is an integer, it lies in the above specified 
range, and Eqs.flph greatly simplify to: 



I, 



= 1/13] b = (3e 



(7) 



where only one of l m ax, (3 and b is a free variable. Choos- 
ing (3 and b such that l max = 20, 30, the result of the 
recursion relation Eq.([2J) is plotted in Fig.l, with nor- 



(3) malisation ^ n(l) = 1. 



Our simulation results are obtained following the 
method proposed by Penna p3, where mutation rate 
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FIG. 2: Mortality rate with 1,, 



30, (note log scale). 



is pre-spccificd as is the maximum string length l S i m (of- 
ten set to 32). However, a Verhulst factor is included 
to regulate the birth reproduction rate and the system 
finds its own stationarity by adjusting the birth rate. 
The string length dependence has been investigated by 
computer simulation, but no universal results were found 
p5[ . Our analysis suggests that the results depend on the 
value of l s i m in relation to l ma x which is set by the mu- 
tation rate, Eq.(Q). If i s j m > l max , no string length effect 
is expected as the final bits of the string do not affect 
the system. Therefore, the most efficient simulation is 
performed when l ma x approaches l S i m from below. 

With the original Penna model, the only mortality at 
age x occurs to individuals with I = x; these have number 
density n(l)/l and mortality rate I. In other words, the 
survival function of an Z-type is a square function with a 
height of I and width I. The normalised mortality at x 
is therefore: 



M{x) 



n(x)/x 



(8) 



where the denominator gives the total number of indi- 
viduals living at the age x, and the numerator gives the 
number of individuals dying at the age x. The resulting 
population mortality exhibits a Gompertz-like behaviour 
up to an intermediate age but then slows down as l max is 
approached j^J, see Fig. 2. But since the mortality at this 
point reaches one, there is no population left to continue 
the incipient plateau and only its onset is observed. 

However, equations (|^-^) make no assumption about 
the survival function of an Z-type and apply equally well 
to, say, a Fermi function [|l7| [l8| : 



exp 



x—l 

w I 



where w is its width in units of the 'Fermi level' I. TV 
ensures /(0,Z) = 1, and I is chosen so that the average 
life span f(x, I) = I. In fact, I = I to within terms of 
order exp[— 1/w] ; the latter are safely ignored in following 
plots, where we choose small w's as examples of modest 
influence of environmental factors on the life expectancy 
of genetically identical individuals. The analysis leading 



FIG. 3: Mortality rate for Fermi survival functions of different 
width w, and l max = 30. 



to the recursion result Eq.(g) stands, and the mortality 
function corresponding to Fig. 2 is plotted in Fig. 3. 

A pronounced 'mortality plateau' is now observed. By 
varying w in a reasonable range of 0.01 — 0.5 and the 
choice of b, which in turn determines n(l), it is possi- 
ble to obtain different shapes of this 'mortality plateau' 
which resemble those observed experimentally in vari- 
ous species || . The different values of w could be easily 
justified as the species-dependent susceptibility to envi- 
ronmental variations as well as other factors. Thus, the 
Penna model, coupled to a Fermi survival function, is suf- 
ficient to account for the plateau observed in mortality 
for different species. 

Now we extend the analysis to the multi-disease cases 
where T > 1. In this scenario, the obvious extension 
of our existing theory is to write n(li, I2, It) which 
gives the number of individuals with deleterious bits at 
positions ?i, Za> —>It on their 'genetic' strings. And a 
similar, albeit more complex, evolution equation can be 
constructed. However, it is immediately evident that the 
final It holds a special position as this is the bit which 
determines the individual's death. In contrast, the posi- 
tions of other deleterious bits are less significant, in fact 
so insignificant as to inspire an ansatz: 

n T (l) = n(h,l 2 , 

for the steady state. Neglecting multiple mutations which 
will be addressed elsewhere ]26j , we find that the relevant 
generalization of Eq.(Q) reduces to: 



= 



dnr(l) 
dt 



= &n T (Z)e-W- T+1 >-^ (9) 



+6T(1 - e H») e -0O-r+i) J2 n T (l'). 

I' =1+1 

Our notation means rii(l) = n(l), defined earlier. The 
ansatz implies that the distribution function n(Zi, I2, I) 
depends merely on the position I and not on any of the 
others. It is possible that there are other solutions not 
contained in the ansatz; without claiming uniqueness we 
note that our simulations show the ansatz leads to a valid 
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Since for every I there are different combinations of hav- 
ing the remaining T — 1 mutations, the correct normal- 
isation is now Y] { Cip_i ut(I) = f, where O r _ 1 is the 
number of combinations of choosing T — 1 out of a total 
of I, and C l T _ x n.r(0 gives the weighted probability of 
finding an i-type. The stability criterion for the multi- 
disease case reads: 

Imax = V/3; b = (3 e 1 -^-^ (10) 
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Mean Lifespan L which recovers Eq.(M) for T = 1. Choosing I = 30 

as before and T — 4, we plot the weighted mutation 
FIG. 4: Lifespan distribution function Ct-iTItCO for T = 4. distribution function in Fig||. The mortality rate with 

the same Fermi survival function as before is presented 
in Fig.| 
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+ + + . vw^wwfflTO.^;: j n conc i us j 011) we have analysed the T = 1 case of 

^^ ++ ++++"^ ^^^^^h . ^ ne p enna model, finding a nontrivial steady state solu- 

tion for an arbitrary i-dependent survival function, and 

^ _ gained by ansatz similar solutions for the general T case. 

Our result suggests the mortality plateau arises from 

a variation in the mortality within genetically identical 

subpopulation. Our analysis also shows the importance 

of the maximum sustainable longevity l max , set by the 

02 — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — 7^ — ^ mutation rate, which determines for example whether 

Age there is any effect of the bit-string length l sim used in 

simulations (there is such an effect only if l S i m < l max ). 

FIG. 5: Mortality rate for T = 4. It is hoped that the analytic solution presented here will 

provide guidance for further simulations. In our own fu- 

, , . t-, a , , , , . ture work we will describe the continuum limit (in I) of 

solution. hq.(H) can be solved to give the recursion: , . ,. . . . ... v ' . 

^ the theory presented here and deal with other, more tech- 

n (I + 1) I + 1 e /3(i-T+i) _ foi nical aspects of the Penna model such as multiple muta- 
ttt — = — : „ : „, — — — — — — ^- tions and the Verhulst factor coupled with noise pfl. 
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